Chapter 6 Does the antibiotic administration change the microbial community?
6.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point %in% c("Acclimation", "Antibiotics") ) label_map <- c(
"richness" = "Species Richness",
"neutral" = "Neutral Diversity",
"phylogenetic" = "Phylogenetic Diversity")
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic")),
time_point = factor(time_point, levels = c("Wild", "Acclimation"))) %>%
filter(time_point %in% c("Wild", "Acclimation")) %>%
ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
geom_boxplot(width = 0.5,alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
geom_jitter(width = 0.2,show.legend = FALSE) +
scale_color_manual(values=c('#008080', "#d57d2c")) +
scale_fill_manual(values=c('#008080', "#d57d2c")) +
facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.text = element_text(size = 14),
axis.title = element_text(size = 14),
strip.text = element_text(size = 18, lineheight = 0.6),
) +
ylab("Alpha diversity")Richness
Modelq0GLMMNB <- glmer.nb(richness ~ time_point*Population+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(5.1645) ( log )
Formula: richness ~ time_point * Population + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
434.8 446.1 -211.4 422.8 43
Scaled residuals:
Min 1Q Median 3Q Max
-1.61979 -0.56608 0.06074 0.51882 2.39923
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.2193 0.4683
Number of obs: 49, groups: individual, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.7459 0.1623 23.073 < 2e-16 ***
time_pointAntibiotics -1.1523 0.1862 -6.190 6.02e-10 ***
PopulationWarm 0.7055 0.2719 2.595 0.00946 **
time_pointAntibiotics:PopulationWarm -0.3586 0.3037 -1.181 0.23769
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) tm_pnA PpltnW
tm_pntAntbt -0.466
PopulatnWrm -0.598 0.278
tm_pntAn:PW 0.292 -0.611 -0.461
Neutral
Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
342.5018 353.3418 -165.2509
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.282619 7.923552
Fixed effects: neutral ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 19.11713 2.078908 25 9.195758 0e+00
time_pointAntibiotics -11.46317 2.832826 20 -4.046551 6e-04
PopulationWarm 24.97345 3.534826 25 7.064972 0e+00
time_pointAntibiotics:PopulationWarm -20.91701 4.793363 20 -4.363745 3e-04
Correlation:
(Intr) tm_pnA PpltnW
time_pointAntibiotics -0.633
PopulationWarm -0.588 0.373
time_pointAntibiotics:PopulationWarm 0.374 -0.591 -0.632
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.8258256 -0.5685215 -0.1867832 0.5823712 2.0651455
Number of Observations: 49
Number of Groups: 27
Phylogenetic
Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
203.6546 214.4946 -95.8273
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.7173285 1.688356
Fixed effects: phylogenetic ~ time_point * Population
Value Std.Error DF t-value p-value
(Intercept) 5.365728 0.4446271 25 12.067930 0.0000
time_pointAntibiotics -0.761455 0.6038652 20 -1.260969 0.2218
PopulationWarm 1.097291 0.7560383 25 1.451370 0.1591
time_pointAntibiotics:PopulationWarm -0.702255 1.0216421 20 -0.687379 0.4997
Correlation:
(Intr) tm_pnA PpltnW
time_pointAntibiotics -0.631
PopulationWarm -0.588 0.371
time_pointAntibiotics:PopulationWarm 0.373 -0.591 -0.629
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.85993684 -0.51436620 -0.05130501 0.48944142 1.80209573
Number of Observations: 49
Number of Groups: 27
6.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Antibiotics")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point %in% c("Acclimation", "Antibiotics"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03530 0.035300 10.073 999 0.002 **
Residuals 47 0.16471 0.003505
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 1.9348458 | 0.1002347 | 6.065418 | 0.001 |
| Population | 1 | 2.1356198 | 0.1106358 | 6.694811 | 0.001 |
| time_point:Population | 1 | 0.8778553 | 0.0454773 | 2.751930 | 0.004 |
| Residual | 45 | 14.3548318 | 0.7436522 | NA | NA |
| Total | 48 | 19.3031527 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.051657 0.051657 9.9224 999 0.005 **
Residuals 47 0.244689 0.005206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 2.0429592 | 0.11061666 | 7.254623 | 0.001 |
| Population | 1 | 2.8764490 | 0.15574623 | 10.214376 | 0.001 |
| time_point:Population | 1 | 0.8770558 | 0.04748846 | 3.114457 | 0.003 |
| Residual | 45 | 12.6723552 | 0.68614864 | NA | NA |
| Total | 48 | 18.4688192 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.67439 0.67439 55.932 999 0.001 ***
Residuals 47 0.56670 0.01206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(phylo ~ time_point*Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
by="terms") %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| time_point | 1 | 1.8783556 | 0.23032747 | 16.358708 | 0.001 |
| Population | 1 | 0.7551942 | 0.09260332 | 6.577030 | 0.001 |
| time_point:Population | 1 | 0.3545686 | 0.04347786 | 3.087959 | 0.023 |
| Residual | 45 | 5.1670341 | 0.63359135 | NA | NA |
| Total | 48 | 8.1551525 | 1.00000000 | NA | NA |
dbRDA
#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(
biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationWarm" = "Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
"Interaction Wam population"
)
CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(values=c('#008080', "#d57d2c")) +
scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(
biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationWarm" = "Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
"Interaction Wam population"
)
CAP_df %>%
group_by(Population, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(values=c('#008080', "#d57d2c")) +
scale_fill_manual(values=c('#00808050', "#d57d2c50"))+
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
rownames_to_column('Tube_code') %>%
left_join(subset_meta, by = 'Tube_code') %>%
column_to_rownames('Tube_code')%>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE))
biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
rownames_to_column("Variable")
biplot_scores$Variable <- recode(
biplot_scores$Variable,
"subset_meta$time_pointAntibiotics" = "Antibiotics",
"subset_meta$PopulationWarm" = "Warm population",
"subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
"Interaction Wam population"
)
CAP_df %>%
group_by(type, time_point) %>%
mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
scale_color_manual(values=c('#008080', "#d57d2c")) +
scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
geom_point(size=2) +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_vline(xintercept = 0, linetype = "dotted") +
geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
theme_classic()6.3 Differences in MAGs and functional capacity
6.3.1 Cold population
6.3.1.1 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold" & time_point == "Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Cold" & time_point == "Antibiotics") %>%
dplyr::select(sample) %>% pull()
subset_meta <- sample_metadata %>%
filter(Population == "Cold" & time_point %in% c("Antibiotics","Acclimation"))
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>%
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 9 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 43
2 p__Bacteroidota 15
3 p__Bacillota 8
4 p__Pseudomonadota 8
5 p__Cyanobacteriota 3
6 p__Verrucomicrobiota 2
7 p__Bacillota_B 1
8 p__Bacillota_C 1
9 p__Fusobacteriota 1
# A tibble: 82 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_1:bin_000006 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__Kineothrix s__
2 AH1_2nd_19:bin_000005 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__Roseburia s__
3 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae g__Akkermansia s__
4 AH1_2nd_2:bin_000003 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__JAAYNV01 s__
5 AH1_2nd_10:bin_000007 p__Bacillota_A c__Clostridia o__Lachnospirales f__ g__ s__
6 AH1_2nd_15:bin_000043 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__
7 AH1_2nd_10:bin_000041 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Rikenellaceae g__Alistipes s__
8 AH1_2nd_8:bin_000007 p__Bacillota c__Bacilli o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma s__
9 AH1_2nd_14:bin_000015 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__
10 AH1_2nd_5:bin_000001 p__Bacillota_C c__Negativicutes o__Selenomonadales f__ g__ s__
# ℹ 72 more rows
Antibiotics
# A tibble: 4 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_7:bin_000001 p__Bacillota_A c__Clostridia o__Clostridiales f__Clostridiaceae g__Sarcina s__
3 AH1_2nd_7:bin_000055 p__Bacillota c__Bacilli o__Mycoplasmatales f__Mycoplasmoidaceae g__Ureaplasma s__
4 AH1_2nd_13:bin_000025 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__
Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>
6.3.1.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold" & time_point %in% c("Acclimation", "Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output_antib_accl_2803 = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output_antib_accl_2803$res %>%
dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
filter(p_time_pointAntibiotics < 0.05) %>%
dplyr::arrange(p_time_pointAntibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_pointAntibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) +
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
MAGs enriched in acclimation on the left side and in antibiotics on the right side
Phyla of the significant MAGs
Acclimation
ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacteroidota 14
2 Bacillota_A 5
3 Bacillota 1
4 Campylobacterota 1
phylum genus species
1 Bacteroidota Bacteroides
2 Campylobacterota Helicobacter_J
3 Bacillota_A
4 Bacteroidota Odoribacter
5 Bacillota
6 Bacteroidota Bacteroides
7 Bacteroidota Bacteroides
8 Bacteroidota Phocaeicola
9 Bacteroidota Bacteroides
10 Bacteroidota Odoribacter
11 Bacteroidota Phocaeicola
12 Bacteroidota Bacteroides
13 Bacteroidota Bacteroides
14 Bacteroidota Parabacteroides
15 Bacteroidota
16 Bacteroidota Parabacteroides
17 Bacillota_A
18 Bacillota_A Lacrimispora
19 Bacteroidota Alistipes
20 Bacillota_A
21 Bacillota_A
Antibiotics
ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics>0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacillota_A 2
2 Bacillota 1
3 Pseudomonadota 1
4 Verrucomicrobiota 1
phylum genus species
1 Bacillota_A MGBC116941
2 Bacillota CAG-345
3 Bacillota_A Intestinimonas
4 Verrucomicrobiota
5 Pseudomonadota
6.3.1.3 Differences in the functional capacity
sample_sub <- sample_metadata %>%
filter(Population == "Cold" & time_point %in% c("Acclimation", "Antibiotics"))
genome_counts_filtered <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample)))
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
genome_counts_filtered_row <- genome_counts_filtered %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)6.3.1.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(time_point) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
time_point MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.335 0.0323
2 Antibiotics 0.247 0.136
MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94311, p-value = 0.09178
Wilcoxon rank sum exact test
data: value by time_point
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(8,9)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(time_point) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c("#003a45",'#008080')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="time_point")6.3.1.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(time_point) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
time_point MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.329 0.0319
2 Antibiotics 0.256 0.126
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95736, p-value = 0.2323
Wilcoxon rank sum exact test
data: value by time_point
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(8,9)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | time_point |
|---|---|---|---|---|---|
| B06 | 0.68059020 | 0.47329940 | Organic anion biosynthesis | 0.20729080 | Acclimation |
| B02 | 0.59930820 | 0.41576970 | Amino acid biosynthesis | 0.18353850 | Acclimation |
| D02 | 0.38530780 | 0.20616790 | Polysaccharide degradation | 0.17913990 | Acclimation |
| S03 | 0.27145162 | 0.09403129 | Spore | 0.17742033 | Acclimation |
| B01 | 0.84159250 | 0.69078390 | Nucleic acid biosynthesis | 0.15080860 | Acclimation |
| B07 | 0.44505530 | 0.30423320 | Vitamin biosynthesis | 0.14082210 | Acclimation |
| B08 | 0.44618700 | 0.32094170 | Aromatic compound biosynthesis | 0.12524530 | Acclimation |
| D09 | 0.25519710 | 0.13446900 | Antibiotic degradation | 0.12072810 | Acclimation |
| B04 | 0.54481600 | 0.42941310 | SCFA biosynthesis | 0.11540290 | Acclimation |
| D03 | 0.29173710 | 0.20687250 | Sugar degradation | 0.08486460 | Acclimation |
| D05 | 0.28853350 | 0.22270070 | Amino acid degradation | 0.06583280 | Acclimation |
| B10 | 0.05960097 | 0.03635236 | Antibiotic biosynthesis | 0.02324861 | Acclimation |
6.3.2 Warm population
6.3.2.1 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Warm" & time_point == "Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Warm" & time_point == "Antibiotics") %>%
dplyr::select(sample) %>% pull()
sample_sub <- sample_metadata %>%
filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics"))
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>%
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 12 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 64
2 p__Bacteroidota 32
3 p__Bacillota 16
4 p__Cyanobacteriota 9
5 p__Pseudomonadota 8
6 p__Bacillota_B 3
7 p__Desulfobacterota 3
8 p__Bacillota_C 2
9 p__Verrucomicrobiota 2
10 p__Actinomycetota 1
11 p__Campylobacterota 1
12 p__Elusimicrobiota 1
# A tibble: 142 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae g__Akkermansia s__
2 LI1_2nd_8:bin_000013 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae g__Akkermansia s__
3 AH1_2nd_19:bin_000005 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__Roseburia s__
4 AH1_2nd_20:bin_000009 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__Kineothrix s__
5 AH1_2nd_18:bin_000011 p__Bacillota c__Bacilli o__Erysipelotrichales f__Coprobacillaceae g__Coprobacillus s__
6 LI1_2nd_2:bin_000001 p__Bacillota_A c__Clostridia o__Oscillospirales f__Ruminococcaceae g__ s__
7 AH1_2nd_6:bin_000035 p__Bacillota c__Bacilli o__RF39 f__UBA660 g__Faecisoma s__
8 LI1_2nd_4:bin_000003 p__Bacillota c__Bacilli o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma s__
9 AH1_2nd_1:bin_000039 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Marinifilaceae g__Odoribacter s__
10 AH1_2nd_2:bin_000018 p__Elusimicrobiota c__Elusimicrobia o__Elusimicrobiales f__Elusimicrobiaceae g__UBA1436 s__
# ℹ 132 more rows
Antibiotics
# A tibble: 7 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_8:bin_000030 p__Actinomycetota c__Actinomycetia o__Mycobacteriales f__Mycobacteriaceae g__Corynebacterium s__
3 LI1_2nd_8:bin_000077 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Tannerellaceae g__Parabacteroides s__Parabacteroides sp003480915
4 LI1_2nd_5:bin_000032 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA1242 g__Caccovivens s__
5 AH1_2nd_18:bin_000013 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Tannerellaceae g__Parabacteroides s__
6 LI1_2nd_5:bin_000069 p__Bacillota c__Bacilli o__RF39 f__UBA660 g__ s__
7 AH1_2nd_20:bin_000061 p__Bacillota c__Bacilli o__Lactobacillales f__Enterococcaceae g__Enterococcus s__
Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 20 × 3
trait p_value p_adjust
<chr> <dbl> <dbl>
1 B0220 0.00618 0.0497
2 B0708 0.00280 0.0497
3 B0709 0.00662 0.0497
4 B1012 0.00618 0.0497
5 D0203 0.00559 0.0497
6 D0205 0.00280 0.0497
7 D0206 0.00280 0.0497
8 D0207 0.00280 0.0497
9 D0210 0.00280 0.0497
10 D0211 0.00662 0.0497
11 D0213 0.00685 0.0497
12 D0305 0.00685 0.0497
13 D0308 0.00685 0.0497
14 D0610 0.00618 0.0497
15 D0706 0.00662 0.0497
16 D0902 0.00618 0.0497
17 D0905 0.00618 0.0497
18 D0908 0.00280 0.0497
19 D0911 0.00618 0.0497
20 S0104 0.00685 0.0497
6.3.2.2 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output_2803 = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output_2803$res %>%
dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
filter(p_time_pointAntibiotics < 0.05) %>%
dplyr::arrange(p_time_pointAntibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_pointAntibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) +
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))
MAGs enriched in acclimation on the left side and in antibiotics on the right side
Phyla of the significant MAGs Acclimation
phylum genus species
1 Bacteroidota Phocaeicola
Antibiotics
ancombc_rand_table_mag%>%
filter(lfc_time_pointAntibiotics>0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacteroidota 2
2 Bacillota 1
3 Bacillota_A 1
4 Chlamydiota 1
5 Pseudomonadota 1
phylum genus species
1 Bacteroidota Bacteroides
2 Bacillota_A MGBC116941
3 Bacteroidota Bacteroides
4 Bacillota Ureaplasma
5 Pseudomonadota
6 Chlamydiota
6.3.2.3 Differences in the functional capacity
sample_sub <- sample_metadata %>%
filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics"))
genome_counts_filtered <- genome_counts_filt %>%
select(one_of(c("genome",sample_sub$sample)))
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
genome_counts_filtered_row <- genome_counts_filtered %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)6.3.2.3.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(time_point) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
time_point MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.355 0.0228
2 Antibiotics 0.253 0.0674
MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.88114, p-value = 0.03317
Wilcoxon rank sum exact test
data: value by time_point
W = 70, p-value = 0.0003291
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(8,9)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(time_point) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
scale_fill_manual(values=c("#003a45",'#008080')) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="time_point")6.3.2.3.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(time_point) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
time_point MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.348 0.0158
2 Antibiotics 0.268 0.0586
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.85791, p-value = 0.01419
Wilcoxon rank sum exact test
data: value by time_point
W = 69, p-value = 0.0005759
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(8,9)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | time_point |
|---|---|---|---|---|---|
| D02 | 0.42716670 | 0.21793910 | Polysaccharide degradation | 0.20922760 | Acclimation |
| B02 | 0.62784110 | 0.47187490 | Amino acid biosynthesis | 0.15596620 | Acclimation |
| D03 | 0.34000130 | 0.18813800 | Sugar degradation | 0.15186330 | Acclimation |
| B08 | 0.44482920 | 0.33165890 | Aromatic compound biosynthesis | 0.11317030 | Acclimation |
| D05 | 0.33300830 | 0.22491860 | Amino acid degradation | 0.10808970 | Acclimation |
| D09 | 0.24738110 | 0.14240830 | Antibiotic degradation | 0.10497280 | Acclimation |
| B06 | 0.70894950 | 0.61231610 | Organic anion biosynthesis | 0.09663340 | Acclimation |
| B09 | 0.02947309 | 0.01568554 | Metallophore biosynthesis | 0.01378755 | Acclimation |